2.1 The shape of the data

## [1] "dimensions: 166" "dimensions: 7"

There are 166 rows with 7 columns each.

The spec is:

## cols(
##   gender = col_character(),
##   age = col_double(),
##   attitude = col_double(),
##   deep = col_double(),
##   stra = col_double(),
##   surf = col_double(),
##   points = col_double()
## )

2.2 Descriptive statistics

Sex distribution

gender seems to be a categorical value, so let’s see the number of rows per gender (sex):

## # A tibble: 2 × 2
##   gender count
##   <chr>  <int>
## 1 F        110
## 2 M         56

110 F, and 56 M, there is a skew towards female students in the dataset. Let’s plot that.

## Loading required package: ggplot2

## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2

fig. 1.1, sex distributions

Age distribution

Let’s plot the age distribution of the students.

lrn2014 %>%
  ggplot(aes(x = age)) +
  stat_count()

fig. 1.2, age distributions

min(lrn2014$age)

## [1] 17

max(lrn2014$age)

## [1] 55

median(lrn2014$age)

## [1] 22

The age range is from 17 to 55, and the median is 22. Visually inspecting the distribution, the mode of the distribution is early twenties, as you would expect, although there is a long tail.

Age-sex distribution

Let’s combine the two columns into a classic population pyramid, or age-sex pyramid.

But that’s not exactly what we want. It turns out a population pyramid is not an out-of-the-box plot we can easily produce, we have to manually calculate bins and do some magic.

lrn2014 %>%
  mutate(age_bin = cut(age, breaks=seq(0,60,5))) %>% # Bin by age
  group_by(gender, age_bin) %>%                      # Group by bin and gender
  summarize(count =n()) %>%                          # Sum over groups
  mutate(count = 
      if_else(gender == 'M', -count, count)) %>%     # Turn one negative
  ggplot(aes(x=count, y = age_bin)) +
  geom_col(aes(fill=gender)) 

## `summarise()` has grouped output by 'gender'. You can override using the
## `.groups` argument.

fig. 1.3a, population pyramid

There are very few male students under 20, I would speculate that this is due to Finnish army conscription, otherwise the distribution seems roughly equal on the female and male sides.

We can of course bin by year instead of 5 years, and we get a higher resolution but more noise.

## `summarise()` has grouped output by 'gender'. You can override using the
## `.groups` argument.

fig. 1.3b, population pyramid #2, fine-grained bins

There is one peculiar decline in the female student participation around ~26-28 which jumps back after thirty. This might be a maternity effect, but this is highly speculative, there’s very few samples in this dataset.

Exam scores

Let’s look at exam scores:

paste("median:", median(lrn2014$points), ", mean:",mean(lrn2014$points), ", standard deviation:", sd(lrn2014$points))

## [1] "median: 23 , mean: 22.7168674698795 , standard deviation: 5.89488364245529"

Let’s look at the full distribution usin geom_density.

lrn2014 %>%
  ggplot(aes(x=points)) +
  geom_density()

fig. 1.4a, exam score distribution

There’s a curious valley in the density at around 12-14 points. Let’s look closer.

lrn2014 %>%
  ggplot(aes(x=points, tickwidth=1)) +
  geom_histogram(boundary=0,binwidth=1) +
  scale_x_continuous(breaks = seq(0, 40, by = 1))

fig. 1.4b, exam score histogram

So no students got 12,13, or 14 points. The jump from 11 to 15 must be behind some barrier. Let’s look at our two demographic variables.

Exploring exam score distributions for different groups

Point means by sex, with age gradient

lrn2014 %>%
  group_by(gender) %>%
  ggplot(aes(x=points)) +
  geom_histogram(boundary=0,binwidth=2, aes(fill=factor(age))) +
  facet_wrap(~gender)

fig. 1.5a, exam scores by sex, with age gradient

I see no clear bias either way, perhaps a slight correlation with score and age within the mode (20-30 years) and then no correlation for higher ages. The female distribution seems most well-behaved, although the gap from 11 points up is much sharper here as well.

Point means by age group, with sex

lrn2014 %>%
  mutate(age_bin = cut(age, breaks=seq(0,60,5))) %>% # Bin by age
  group_by(gender, age_bin) %>%
  ggplot(aes(x=points)) +
  geom_histogram(boundary=0,binwidth=2, aes(fill=gender)) +
  facet_wrap(~age_bin)

fig. 1.5b, exam scores by age group, with labeled sex

Point means by age group, box plot

lrn2014 %>%
  mutate(age_bin = cut(age, breaks=seq(0,60,5))) %>% # Bin by age
  ggplot(aes(x=age_bin, y = points)) +
  stat_boxplot(geom = "errorbar", width = 0.25) +
      geom_boxplot()

fig. 1.6a, exam score distributions by age group, box plot sequence

lrn2014 %>%
  mutate(age_bin = cut(age, breaks=seq(0,60,5))) %>% # Bin by age
  ggplot(aes(x=age_bin, y = points)) +
  stat_boxplot(geom = "errorbar", width = 0.25) +
      geom_boxplot()+
  facet_wrap(~gender)

fig. 1.6b, exam score distributions by age-sex group, box plot sequence

I see no correlation between age and exam score in these plots.

2.3 Correlation matrix

Finally, let us look at all the survey question scores together with the already explored variables in a correlation matrix.

lrn2014 %>%
  select(gender, age, surf, stra, deep, attitude, points) %>%
  ggpairs(aes(fill=gender, color=gender),lower=list(combo=wrap("facethist",binwidth=0.5))) # Have to add the lower arg so ggpairs doesn't complain

fig. 1.7, correlation matrix of all variables

There seem to be negative correlations between: - surf and deep (but seemingly only strongly for male students) - attitude and deep (weak, but stronger for male students again) - stra and surf (weak)

And a possitive correlation between points and attitude! This is the strongest linear relationship in the data. And we can verify that age does not seem to have an effect at all.

There also seems to be a relationship between attitude and gender, but no relationship between attitude and points.

3.1 Simple regression as baseline

Let’s start with a simple regression model of points ~ attitude, which will be our baseline.

attitude_model <- lrn2014 %>%
                  lm(points ~ attitude, data = .)
summary(attitude_model)

## 
## Call:
## lm(formula = points ~ attitude, data = .)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -16.9763  -3.2119   0.4339   4.1534  10.6645 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  11.6372     1.8303   6.358 1.95e-09 ***
## attitude      3.5255     0.5674   6.214 4.12e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.32 on 164 degrees of freedom
## Multiple R-squared:  0.1906, Adjusted R-squared:  0.1856 
## F-statistic: 38.61 on 1 and 164 DF,  p-value: 4.119e-09

There is clearly a statistically significant relationship between attitude and points. But R-squared is only around 18.5%, so there is a lot of variance not explained by the model.

3.2 Multiple regression, 3 variables

three_var_model <- lrn2014 %>%
                  lm(points ~ attitude + stra + surf, data = .)
summary(three_var_model)

## 
## Call:
## lm(formula = points ~ attitude + stra + surf, data = .)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -17.1550  -3.4346   0.5156   3.6401  10.8952 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  11.0171     3.6837   2.991  0.00322 ** 
## attitude      3.3952     0.5741   5.913 1.93e-08 ***
## stra          0.8531     0.5416   1.575  0.11716    
## surf         -0.5861     0.8014  -0.731  0.46563    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.296 on 162 degrees of freedom
## Multiple R-squared:  0.2074, Adjusted R-squared:  0.1927 
## F-statistic: 14.13 on 3 and 162 DF,  p-value: 3.156e-08

The adjusted R-squared is higher, which means our model is capturing more of the underlying interactions than before, although still below 20%.

It seems that the relationship between points and surf is not statistically significant.

Let’s drop surf and keep stra, and try again.

3.3 Multiple regression, 2 variables

two_var_model <- lrn2014 %>%
                  lm(points ~ attitude + stra, data = .)
summary(two_var_model)

## 
## Call:
## lm(formula = points ~ attitude + stra, data = .)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -17.6436  -3.3113   0.5575   3.7928  10.9295 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   8.9729     2.3959   3.745  0.00025 ***
## attitude      3.4658     0.5652   6.132 6.31e-09 ***
## stra          0.9137     0.5345   1.709  0.08927 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.289 on 163 degrees of freedom
## Multiple R-squared:  0.2048, Adjusted R-squared:  0.1951 
## F-statistic: 20.99 on 2 and 163 DF,  p-value: 7.734e-09

Right, this is our best fit yet, based on the adjusted R-squared. Depending on what our choice of significance level would be in a hypothesis test, the interaction with stra would be ignored. At a standard level of a = 0.05, we wouldn’t reject the null hypothesis that stra has a linear effect on points.

Let’s test another model, with nothing but stra.

stra_model <- lrn2014 %>%
                  lm(points ~ stra, data = .)
summary(stra_model)

## 
## Call:
## lm(formula = points ~ stra, data = .)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -16.5581  -3.8198   0.1042   4.3024  10.1394 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   19.233      1.897  10.141   <2e-16 ***
## stra           1.116      0.590   1.892   0.0603 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.849 on 164 degrees of freedom
## Multiple R-squared:  0.02135,    Adjusted R-squared:  0.01538 
## F-statistic: 3.578 on 1 and 164 DF,  p-value: 0.06031

3.4 Interpretation of the statistical test related to the model

The statistical tests that are reported in the summary are:

• the t-test
• the F-test

4.1 Summary

Let’s rerun that summary:

summary(attitude_model)

## 
## Call:
## lm(formula = points ~ attitude, data = .)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -16.9763  -3.2119   0.4339   4.1534  10.6645 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  11.6372     1.8303   6.358 1.95e-09 ***
## attitude      3.5255     0.5674   6.214 4.12e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.32 on 164 degrees of freedom
## Multiple R-squared:  0.1906, Adjusted R-squared:  0.1856 
## F-statistic: 38.61 on 1 and 164 DF,  p-value: 4.119e-09

The fitted regression coefficients are: - intercept: 11.6372 - attitude: 3.5255

Which means the that the model predicts the conditional mean of the exam scores, given an attitude value as:

points = 11.6372 + 3.5255*attitude

If we multiply the attitude coefficient with the range of the attitude in the population, we can get an idea of how the model assigns expected exam scores based on a students attitude.

am <- mean(lrn2014$attitude)
as <- sd(lrn2014$attitude)
print("Range of predictor term within a sd:")

## [1] "Range of predictor term within a sd:"

3.5255*c(am-as, am, am+as)

## [1]  8.506549 11.079839 13.653130

print("Range of ŷ_i within a sd:")

## [1] "Range of ŷ_i within a sd:"

11.6732 + 3.5255*c(am-as, am, am+as)

## [1] 20.17975 22.75304 25.32633

print("Range of ŷ_i within two sd's:")

## [1] "Range of ŷ_i within two sd's:"

11.6732 + 3.5255*c(am-2*as, am, am+2*as)

## [1] 17.60646 22.75304 27.89962

spec(lrn2014)

## cols(
##   gender = col_character(),
##   age = col_double(),
##   attitude = col_double(),
##   deep = col_double(),
##   stra = col_double(),
##   surf = col_double(),
##   points = col_double()
## )

So, assuming attitude is Gaussian and that the sample stddev is a good estimate, the model assigns: an exam score: - between 20.17975 and 25.32633 to a majority of students (about 68%, one stddev in a Gaussian captures 34.1% of the population) - between 17.60646 and 27.89962 to a super-majority (95%) of students

If we look back at the exam score distribution in figure 1.4, this does capture the mode of the distribution.

3.2 Multiple R-squared

The multiple R-squared value is 0.1906, the standard way to express this is that “19% of the variation in exam scores is explained by the attitude variable” (see MABS4IODS, 3.2.1).

I would interpret this to mean that, using the linear model, given attitude we estimate the expectation of the standard error (squared) of this prediction to be roughly 80% of what it would be when simply using the sample mean as a predictor. The estimation assumes the sample is representative, because we’re using residuals to get this number.

I’m not quite sure if this more elaborate interpretation is exact, but it’s what I was able to piece together from sources online, mostly wikipedia (https://en.wikipedia.org/wiki/Coefficient_of_determination).

5.1 Assumptions

The assumptions of the model are that: 1. the constant variance assumption (homoscedasticity) 2. the normality assumption: 3. there’s a linear relationship between the predictor and the response

attitude_model_bptest <- bptest(attitude_model)
attitude_model_bptest

## 
##  studentized Breusch-Pagan test
## 
## data:  attitude_model
## BP = 0.0039163, df = 1, p-value = 0.9501

shapiro.test(attitude_model[["residuals"]])

## 
##  Shapiro-Wilk normality test
## 
## data:  attitude_model[["residuals"]]
## W = 0.97394, p-value = 0.003195

5.1 Odds ratio

coeffs <- coefficients(model)
odds_ratio <- coeffs %>% exp
odds_ratio

## (Intercept)          G3    absences 
##   0.6788751   0.9267616   1.0878762

The odds ratios for G3 and absences are roughly 0.926 and 1.088 For G3, the OR is less than one because the correlation between the variables is negative.

The way to interpret these is that: - for each decrease in final grade average, the odds increase by roughly 8.0% that a student has high alcohol use (1/0.9267 ~= 1.08) - for each increase in absences per course, the odds increase by roughly 8.7% that a student has high alcohol use

But that’s just the average effect the model has fit, let’s look at confidence intervals in the odds ratio

5.2 Confidence intervals

ci <- confint(model) %>% exp

## Waiting for profiling to be done...

ci

##                 2.5 %    97.5 %
## (Intercept) 0.2857593 1.5832545
## G3          0.8639686 0.9935498
## absences    1.0419618 1.1407755

The 95% confidence intervals are both on one side of unity, so we can say with 95% certainty that there is an effect for both variables, and the effect increases the odds by: - a factor in the range of [1.007, 1.15] for each decrease in final grade average (again,inverting the odds ratio because the effect is negative) - a factor in the range of [1.04, 1.14] for each increase in absences per course

5.3 Interpretation

There is an effect, although the final grade average effect goes damn near 1 at the low end of the confidence interval, the unit increase in odds is only 0.7%! For absences, the confidence interval is a little tighter, but it still seems a little wide to me for practical use. We would need more samples in order to get a tighter interval.

The two hypotheses that survived the initial exploration both match the outcomes of the logistic regression, and now we can quantify the quality of an estimate of the effect.

Caveat

We haven’t done a split into a model fit dataset and a validation set, so this confusion matrix is of limited value.

Dataset structure

str(Boston)

## 'data.frame':    506 obs. of  14 variables:
##  $ crim   : num  0.00632 0.02731 0.02729 0.03237 0.06905 ...
##  $ zn     : num  18 0 0 0 0 0 12.5 12.5 12.5 12.5 ...
##  $ indus  : num  2.31 7.07 7.07 2.18 2.18 2.18 7.87 7.87 7.87 7.87 ...
##  $ chas   : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ nox    : num  0.538 0.469 0.469 0.458 0.458 0.458 0.524 0.524 0.524 0.524 ...
##  $ rm     : num  6.58 6.42 7.18 7 7.15 ...
##  $ age    : num  65.2 78.9 61.1 45.8 54.2 58.7 66.6 96.1 100 85.9 ...
##  $ dis    : num  4.09 4.97 4.97 6.06 6.06 ...
##  $ rad    : int  1 2 2 3 3 3 5 5 5 5 ...
##  $ tax    : num  296 242 242 222 222 222 311 311 311 311 ...
##  $ ptratio: num  15.3 17.8 17.8 18.7 18.7 18.7 15.2 15.2 15.2 15.2 ...
##  $ black  : num  397 397 393 395 397 ...
##  $ lstat  : num  4.98 9.14 4.03 2.94 5.33 ...
##  $ medv   : num  24 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 ...

There are 506 rows with 14 columns. The dataset seems to originally have been put together to analyze housing values across the suburbs. From the paper cited on the stat.ethz.ch site: “Harrison, D. and Rubinfeld, D.L. (1978) Hedonic housing prices and the demand for clean air”, we can find that a row represents a “SMSA census tract”, so, a census area in Boston containing some number of housing units.

The columns contain social statistics related to these census areas (e.g. crim = crime rate, ptratio = pupil-teacher ratio), data about the housing units in the area (rm = avg # of rooms per unit, medv = median housing unit value, age = prop. houses built before 1940), and data about the location of the area (e.g. dis = weighted mean of distances to employment centers, chas = 1 if by Charles River, 0 if not by Charles River).

Some of the columns are a little counter-intuitive or difficult to interpret. E.g., the column age is the proportion of houses built before 1940, and the column lstat is the proportion of the population that is lower status. From the Harrison & Rubinfield paper, lower status means: “1/2 * (proportion of adults without some hig h school education and proportion of male workers classified as laborers)”.

Ok, before we move forward, we did see a small issue here, let’s change chas from an integer to a boolean.

library(dplyr)
Boston_explore <- Boston %>% mutate(chas = chas == 1)
str(Boston_explore)

## 'data.frame':    506 obs. of  14 variables:
##  $ crim   : num  0.00632 0.02731 0.02729 0.03237 0.06905 ...
##  $ zn     : num  18 0 0 0 0 0 12.5 12.5 12.5 12.5 ...
##  $ indus  : num  2.31 7.07 7.07 2.18 2.18 2.18 7.87 7.87 7.87 7.87 ...
##  $ chas   : logi  FALSE FALSE FALSE FALSE FALSE FALSE ...
##  $ nox    : num  0.538 0.469 0.469 0.458 0.458 0.458 0.524 0.524 0.524 0.524 ...
##  $ rm     : num  6.58 6.42 7.18 7 7.15 ...
##  $ age    : num  65.2 78.9 61.1 45.8 54.2 58.7 66.6 96.1 100 85.9 ...
##  $ dis    : num  4.09 4.97 4.97 6.06 6.06 ...
##  $ rad    : int  1 2 2 3 3 3 5 5 5 5 ...
##  $ tax    : num  296 242 242 222 222 222 311 311 311 311 ...
##  $ ptratio: num  15.3 17.8 17.8 18.7 18.7 18.7 15.2 15.2 15.2 15.2 ...
##  $ black  : num  397 397 393 395 397 ...
##  $ lstat  : num  4.98 9.14 4.03 2.94 5.33 ...
##  $ medv   : num  24 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 ...

Summary

Show a graphical overview of the data and show summaries of the variables in the data. Describe and interpret the outputs, commenting on the distributions of the variables and the relationships between them. (0-2 points)

summary(Boston_explore)

##       crim                zn             indus          chas        
##  Min.   : 0.00632   Min.   :  0.00   Min.   : 0.46   Mode :logical  
##  1st Qu.: 0.08205   1st Qu.:  0.00   1st Qu.: 5.19   FALSE:471      
##  Median : 0.25651   Median :  0.00   Median : 9.69   TRUE :35       
##  Mean   : 3.61352   Mean   : 11.36   Mean   :11.14                  
##  3rd Qu.: 3.67708   3rd Qu.: 12.50   3rd Qu.:18.10                  
##  Max.   :88.97620   Max.   :100.00   Max.   :27.74                  
##       nox               rm             age              dis        
##  Min.   :0.3850   Min.   :3.561   Min.   :  2.90   Min.   : 1.130  
##  1st Qu.:0.4490   1st Qu.:5.886   1st Qu.: 45.02   1st Qu.: 2.100  
##  Median :0.5380   Median :6.208   Median : 77.50   Median : 3.207  
##  Mean   :0.5547   Mean   :6.285   Mean   : 68.57   Mean   : 3.795  
##  3rd Qu.:0.6240   3rd Qu.:6.623   3rd Qu.: 94.08   3rd Qu.: 5.188  
##  Max.   :0.8710   Max.   :8.780   Max.   :100.00   Max.   :12.127  
##       rad              tax           ptratio          black       
##  Min.   : 1.000   Min.   :187.0   Min.   :12.60   Min.   :  0.32  
##  1st Qu.: 4.000   1st Qu.:279.0   1st Qu.:17.40   1st Qu.:375.38  
##  Median : 5.000   Median :330.0   Median :19.05   Median :391.44  
##  Mean   : 9.549   Mean   :408.2   Mean   :18.46   Mean   :356.67  
##  3rd Qu.:24.000   3rd Qu.:666.0   3rd Qu.:20.20   3rd Qu.:396.23  
##  Max.   :24.000   Max.   :711.0   Max.   :22.00   Max.   :396.90  
##      lstat            medv      
##  Min.   : 1.73   Min.   : 5.00  
##  1st Qu.: 6.95   1st Qu.:17.02  
##  Median :11.36   Median :21.20  
##  Mean   :12.65   Mean   :22.53  
##  3rd Qu.:16.95   3rd Qu.:25.00  
##  Max.   :37.97   Max.   :50.00

Looking at this summary, all columns are definitely not created equal. I am mostly looking at the difference between the mean and median, which — if they differ by much — can indicate that a variable is not normally distributed. Some of the columns are close to normal distributions, but e.g. zn has a median of 0 and a mean of 11.36. Other highly skewed columns are rad, crim, tax, chas, and black.

With ggpairs

library(ggplot2)
library(GGally)
Boston_explore %>%
  ggpairs(lower=list(combo=wrap("facethist",binwidth=0.5))) # Have to add the lower arg so ggpairs doesn't complain

fig. 3.1, Correlation matrix

Ok, lot’s to unpack here, let’s start with the a visual check of each variable’s distribution (the diagonal). Almost none of the columns look normally distributed, with perhaps the exception of rm, the number of rooms.

There are lots of interesting correlations, just looking at the scatter plots, the three values rm, medv, and lstat seem to have strikingly strong relationships with each other with medv, which makes sense to me.

The rad variable, which is an “index of accessibility to radial highways, is clearly a bimodal, or one could even say a split distribution. A subset of areas have a much higher index than the others, and in the scatter plots, this clearly visible line of that higher-index population seems to consistently cover different ranges of the other variable than the lower-index population. The effect is most clearly noticeable in the crim, tax, nox, dis and black scatter plots.

dis and nox also have a strikingly logarithmic-looking relationship.

In general, nearly every variable seems to be correlated with every other variable, excepting the chas (area is by the Charles river) column.

With corrplot

library(corrplot)

## corrplot 0.92 loaded

corrplot(cor(Boston_explore), method="circle")

fig. 3.2, Correlation matrix with corrplot

Using corrplot, we lose some information, but get a better overview of where the correlations are strongest.

We see strong correlations (large balls) between: - dis and (zn, indus, nox, and age) - tax and rad, and this is a very strong correltaion, they seem to capture much of the same variation within them - tax and (crim, indus, and nox) - ditto for rad - lstat and medv

Create a categorical crime rate column

# Create quartile class
Boston_scaled$crim <- as.numeric(Boston_scaled$crim)
Boston_scaled$crime <- cut(Boston_scaled$crim, breaks = quantile(Boston_scaled$crim), include.lowest = TRUE, label=c("low", "med_low", "med_high", "high"))

# Drop crim
Boston_scaled <- dplyr::select(Boston_scaled, -crim)

We’ve split the crime rate column into a categorical variable defining in which quartile of the crime rate distribution the sensus area is in.

Create training and test datasets

set.seed(179693716) 
n <- nrow(Boston_scaled)
split <- sample(n,  size = n * 0.8)

train <- Boston_scaled[split,]
test <- Boston_scaled[-split,]

Now we’ve split the dataset into two: 80% of rows are in the training set, 20% are in the test set.

Biplot

arrows <- function(x, scale = 1, arrow_heads = 0.1, color = "red", tex = 0.75, choices = c(1,2)){
  heads <- coef(x)
  graphics::arrows(x0 = 0, y0 = 0, 
         x1 = scale * heads[,choices[1]], 
         y1 = scale * heads[,choices[2]], col=color, length = arrow_heads)
  text(scale * heads[,choices], labels = row.names(heads), 
       cex = tex, col=color, pos=3)
}

classes <- as.numeric(train$crime)

plot(lda.fit, dimen = 2, col=classes, pch=classes)
arrows(lda.fit, scale = 1.5, color = "#ee8855")

fig. 5.1 LDA biplot

We see that out of the two first linear discriminants, LD1 nearly perfectly separates the data into two clusters: those with high crime rate, and those with other values. rad has the clearly highest coefficient in LD1, which can be seen both from the biplot and the LDA summary.

LD2 seems to find another axis within the data that explains a smaller effect. The largest coefficients in LD2 belong to nox, medv, and zn.

6.Extra: let’s try with only rad

lda.radfit <- lda(crime ~ rad, data = train)
lda.radpred <- predict(lda.radfit, newdata = test)
radtab <- table(correct = facit, predicted = lda.radpred$class)
image(radtab)

fig. 6.2 Confusion matrix with a single variable LDA fit

Using only rad gives us a model that seems to have exactly the same predictive power in the high vs not high case, but loses information in the lower quartiles. This fits with our earlier analysis of how LD1 was mostly rad and was able to carve out most of the high crime rate areas from the rest.

7.Extra Let’s try with k=3

Why not try with k=3 for good measure? Maybe we can find additional structure in the data.

set.seed(124293)
k3m <- kmeans(Boston_kmeans, centers = 3)
summary(k3m)

##              Length Class  Mode   
## cluster      506    -none- numeric
## centers       42    -none- numeric
## totss          1    -none- numeric
## withinss       3    -none- numeric
## tot.withinss   1    -none- numeric
## betweenss      1    -none- numeric
## size           3    -none- numeric
## iter           1    -none- numeric
## ifault         1    -none- numeric

k3m$centers

##         crim         zn      indus        chas        nox         rm        age
## 1 -0.4076669  1.1526549 -0.9877755 -0.10115080 -0.9634859  0.7739125 -1.1241828
## 2 -0.3688324 -0.3935457 -0.1369208  0.07398993 -0.1662087 -0.1700456  0.1673019
## 3  0.8942488 -0.4872402  1.0913679 -0.01330932  1.1109351 -0.4609873  0.7828949
##           dis        rad        tax     ptratio      black       lstat
## 1  1.05650031 -0.5965522 -0.6837494 -0.62478941  0.3607235 -0.90904433
## 2 -0.07766431 -0.5799077 -0.5409630 -0.04596655  0.2680397 -0.05818052
## 3 -0.84882600  1.3656860  1.3895093  0.63256391 -0.7083974  0.90799414
##          medv
## 1  0.84137443
## 2 -0.04811607
## 3 -0.69550394

ggpairs(Boston_kmeans[c("crim", "rad", "tax", "black", "age", "medv", "dis", "zn", "rm", "lstat")], aes(color=factor(k3m$cluster), fill=factor(k3m$cluster)), lower=list(combo=wrap("facethist",binwidth=0.5)))

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

fig. 7.1 k-means pairs

Here we see that the blue cluster in this plot is roughly the same as the blue cluster from the k=2 clustering. The k=2 red cluster has here been split into red and green.

The differences in the red and green clusters seem to be: - the red cluster has higher values of zn more big plots - the red cluster has a smaller proportion of older buildings - the red cluster consists of exclusively black neighbourhoods, while the green cluster has some spread - the green cluster is between the red and blue clusters when it comes to proportion of laborers and uneducated adults

Maybe the red cluster matches better with those black neighbourhoods built more recently, which the 1949 Housing Act and the Federal Housing Authority regulations apply to? I don’t know for sure, more analysis would be required.

0. Setup, read the data

library(readr)
library(tibble)
library(GGally)
library(corrplot)
library(dplyr)
library(FactoMineR)
human0 <- read_csv("https://raw.githubusercontent.com/KimmoVehkalahti/Helsinki-Open-Data-Science/master/datasets/human2.csv", show_col_types = FALSE)

Rownames

human <- column_to_rownames(human0, "Country")

Summary, statistics

Let’s look at the data, starting with some descriptive statistics:

summary(human)

##     Edu2.FM          Labo.FM          Life.Exp        Edu.Exp     
##  Min.   :0.1717   Min.   :0.1857   Min.   :49.00   Min.   : 5.40  
##  1st Qu.:0.7264   1st Qu.:0.5984   1st Qu.:66.30   1st Qu.:11.25  
##  Median :0.9375   Median :0.7535   Median :74.20   Median :13.50  
##  Mean   :0.8529   Mean   :0.7074   Mean   :71.65   Mean   :13.18  
##  3rd Qu.:0.9968   3rd Qu.:0.8535   3rd Qu.:77.25   3rd Qu.:15.20  
##  Max.   :1.4967   Max.   :1.0380   Max.   :83.50   Max.   :20.20  
##       GNI            Mat.Mor         Ado.Birth         Parli.F     
##  Min.   :   581   Min.   :   1.0   Min.   :  0.60   Min.   : 0.00  
##  1st Qu.:  4198   1st Qu.:  11.5   1st Qu.: 12.65   1st Qu.:12.40  
##  Median : 12040   Median :  49.0   Median : 33.60   Median :19.30  
##  Mean   : 17628   Mean   : 149.1   Mean   : 47.16   Mean   :20.91  
##  3rd Qu.: 24512   3rd Qu.: 190.0   3rd Qu.: 71.95   3rd Qu.:27.95  
##  Max.   :123124   Max.   :1100.0   Max.   :204.80   Max.   :57.50

Based on the difference between mean and median, and the ranges of the variables, Edu.Exp and Parli.FM seem the most normally distributed. Edu2.FM, Labo.FM and Life.Exp are more skewed or have long tails (bdifference between mean and median, and the ranges of the variables). GNI and Mat.Mor are closer to log-normal (or a related distribution), and maybe Ado.Birth as well.

Summary, graphical

corrplot(cor(human), method="circle",type = "upper", tl.pos = "d")

fig. 1.1, correlation matrix of variables

There are definitely strong correlations in the data. We can identify four different sets of variables from this matrix.

First, there is one set of very strongly correlated variables: Mat.Mor, Ado.Birth, Life.Exp, and Edu.Exp. These will likely be very strongly represented in the first principal component of a PCA.

Second, these four variables also correlate very strongly with GNI and Edu2.FM, more strongly than that pair of variables correlate between themselves. Since they correlate less amongst themselves, there seems to be more degrees of freedom here.

Then there is a third group, Labo.FM and Parli.F. They seem much less correlated compared to the other set of variables, however, they are most strongly correlated with each other. Again, seemingly more degrees of freedom in this group.

ggpairs(human, progress = FALSE)

fig. 1.2, distributions and scatterplots of variables

The correlations become much more clear when looking at the scatterplots. Some are clearly linear, like Edu.Exp and Life.Exp. It’s even clearer that GNI might need to be log-transformed, based on the distribution and the scatter plot shapes. Possibly Mat.Mor and Ado.Birth too.

The other distributions look like skewed normal-like distributions, all of them have one big main mode, although there are some interesting concentrations in certain parts of the distributions.

E.g., in Life.Exp, there is a clear plateau from 60 to roughly 65 or so. This may be related to quality of life before and after retirement.

Another example is this curious smaller mode in the Edu2.FM around 0.5, and a subsequent drop after that. This suggests that there are two overlaid populations, one where the mode is a little under 1, another where the mode is around 0.5. I don’t have a hypothesis for what could be causing this difference in behaviors in the two populations of countries, but it’s a very interesting feature of the distribution. Unfortunately we’ve already thrown away the two variables that this variable is based on, which may have given us some clues to what this effect is.

Let’s do some transforms, just for fun

human.log <- human %>% mutate(GNI = log(GNI)) %>% mutate(Mat.Mor = log(Mat.Mor)) %>% mutate(Ado.Birth = log(Ado.Birth))
ggpairs(human.log, progress = FALSE)

fig. 1.3, distributions and scatterplots of variables after log transform

The scatter plots now seem to form much clearer trend lines, which seems to indicate this was a good idea. Let’s recheck the correlation matrix, just to check that the scale of the data didn’t mix things up too much.

corrplot(cor(human.log), method="circle",type = "upper", tl.pos = "d")

fig. 1.4, correlation matrix using log transformed dataset

Well, this does seem to have changed a lot of things. The correlations are now much more even among five variables: Life.Exp, Edu.Exp, GNI, Mat.Mor, and Ado.Birth. The correlation with Edu2.FM is a little weaker, and the Parli.F and Labo.FM variables are again least correlated with the others.

3.1 Interpretation

5.1 MCA

There are a lot of variables here, so let’s choose a sensible subset of them to look at. Let’s choose the following variables, which I’m interpreting free-hand here, since I can’t find much metadata:

• Tea: what kind of tea out of three types the respondent prefers (green, black, Earl Grey)
• price: the price level of the tea that the respondent prefers
• How: whether the respondent takes tea as is, with lemon, with milk, or in some other way
• sex: the sex of the respondent
• SPC: some kind of general social group for the respondent (student, employee (white collar?), middle (management?), senior)
• age_Q: the age group of the respondent
• frequency: how often the respondent drinks tea

of these, Tea, price, and How are “active” variables, and the rest are “supplementary” variables.

As a sidenote, this dataset is very badly documented.

tea_filtered <- tea %>% dplyr::select(one_of("Tea", "price", "How", "sex", "SPC", "age_Q", "frequency"))
summary(tea_filtered)

##         Tea                  price        How      sex               SPC    
##  black    : 74   p_branded      : 95   alone:195   F:178   employee    :59  
##  Earl Grey:193   p_cheap        :  7   lemon: 33   M:122   middle      :40  
##  green    : 33   p_private label: 21   milk : 63           non-worker  :64  
##                  p_unknown      : 12   other:  9           other worker:20  
##                  p_upscale      : 53                       senior      :35  
##                  p_variable     :112                       student     :70  
##                                                            workman     :12  
##    age_Q          frequency  
##  +60  :38   +2/day     :127  
##  15-24:92   1 to 2/week: 44  
##  25-34:69   1/day      : 95  
##  35-44:40   3 to 6/week: 34  
##  45-59:61                    
##                              
## 

Let’s look at a summary of the MCA of the dataset based on these specs above.

mca=MCA(tea_filtered,quali.sup=4:7 ,graph=F)
summary(mca)

## 
## Call:
## MCA(X = tea_filtered, quali.sup = 4:7, graph = F) 
## 
## 
## Eigenvalues
##                        Dim.1   Dim.2   Dim.3   Dim.4   Dim.5   Dim.6   Dim.7
## Variance               0.423   0.399   0.383   0.355   0.335   0.330   0.316
## % of var.             12.688  11.967  11.497  10.653  10.043   9.907   9.484
## Cumulative % of var.  12.688  24.655  36.152  46.805  56.848  66.755  76.239
##                        Dim.8   Dim.9  Dim.10
## Variance               0.283   0.271   0.238
## % of var.              8.490   8.142   7.129
## Cumulative % of var.  84.729  92.871 100.000
## 
## Individuals (the 10 first)
##                    Dim.1    ctr   cos2    Dim.2    ctr   cos2    Dim.3    ctr
## 1               |  0.720  0.409  0.056 |  1.341  1.503  0.196 |  0.789  0.541
## 2               |  0.782  0.481  0.216 |  0.593  0.294  0.124 |  0.071  0.004
## 3               | -0.462  0.168  0.231 |  0.110  0.010  0.013 | -0.707  0.435
## 4               | -0.462  0.168  0.231 |  0.110  0.010  0.013 | -0.707  0.435
## 5               | -0.462  0.168  0.231 |  0.110  0.010  0.013 | -0.707  0.435
## 6               | -0.357  0.100  0.027 | -0.078  0.005  0.001 |  0.344  0.103
## 7               | -0.462  0.168  0.231 |  0.110  0.010  0.013 | -0.707  0.435
## 8               |  0.782  0.481  0.216 |  0.593  0.294  0.124 |  0.071  0.004
## 9               |  0.497  0.195  0.083 |  0.222  0.041  0.016 |  0.544  0.257
## 10              |  1.104  0.961  0.443 | -0.678  0.384  0.167 |  0.138  0.016
##                   cos2  
## 1                0.068 |
## 2                0.002 |
## 3                0.542 |
## 4                0.542 |
## 5                0.542 |
## 6                0.025 |
## 7                0.542 |
## 8                0.002 |
## 9                0.099 |
## 10               0.007 |
## 
## Categories (the 10 first)
##                     Dim.1     ctr    cos2  v.test     Dim.2     ctr    cos2
## black           |   1.366  36.278   0.611  13.516 |  -0.083   0.143   0.002
## Earl Grey       |  -0.440   9.794   0.348 -10.207 |   0.312   5.225   0.175
## green           |  -0.493   2.106   0.030  -2.996 |  -1.636  24.612   0.331
## p_branded       |  -0.497   6.176   0.115  -5.856 |  -0.218   1.257   0.022
## p_cheap         |   1.116   2.289   0.030   2.982 |   1.313   3.361   0.041
## p_private label |   0.018   0.002   0.000   0.084 |  -0.132   0.102   0.001
## p_unknown       |   0.314   0.310   0.004   1.108 |   2.953  29.139   0.363
## p_upscale       |   1.063  15.724   0.242   8.512 |  -0.874  11.272   0.164
## p_variable      |  -0.188   1.036   0.021  -2.504 |   0.225   1.575   0.030
## alone           |  -0.275   3.870   0.140  -6.477 |  -0.328   5.831   0.199
##                  v.test     Dim.3     ctr    cos2  v.test  
## black            -0.825 |   0.099   0.209   0.003   0.977 |
## Earl Grey         7.240 |  -0.247   3.420   0.110  -5.741 |
## green            -9.947 |   1.224  14.341   0.185   7.443 |
## p_branded        -2.566 |   0.459   5.803   0.098   5.403 |
## p_cheap           3.509 |   0.315   0.201   0.002   0.841 |
## p_private label  -0.626 |   1.039   6.569   0.081   4.928 |
## p_unknown        10.421 |   1.519   8.030   0.096   5.362 |
## p_upscale        -6.999 |   0.310   1.476   0.021   2.483 |
## p_variable        2.999 |  -0.913  27.080   0.497 -12.188 |
## alone            -7.721 |  -0.153   1.330   0.044  -3.614 |
## 
## Categorical variables (eta2)
##                   Dim.1 Dim.2 Dim.3  
## Tea             | 0.611 0.359 0.207 |
## price           | 0.324 0.559 0.565 |
## How             | 0.334 0.279 0.378 |
## 
## Supplementary categories (the 10 first)
##                    Dim.1   cos2 v.test    Dim.2   cos2 v.test    Dim.3   cos2
## F               | -0.099  0.014 -2.071 | -0.001  0.000 -0.018 | -0.092  0.012
## M               |  0.145  0.014  2.071 |  0.001  0.000  0.018 |  0.135  0.012
## employee        | -0.153  0.006 -1.307 | -0.104  0.003 -0.891 |  0.027  0.000
## middle          |  0.179  0.005  1.213 |  0.095  0.001  0.645 |  0.033  0.000
## non-worker      |  0.281  0.021  2.531 | -0.214  0.012 -1.927 | -0.010  0.000
## other worker    |  0.221  0.003  1.021 |  0.174  0.002  0.804 |  0.108  0.001
## senior          |  0.226  0.007  1.419 |  0.086  0.001  0.543 |  0.009  0.000
## student         | -0.350  0.037 -3.343 |  0.108  0.004  1.027 | -0.090  0.002
## workman         | -0.327  0.004 -1.154 |  0.167  0.001  0.588 |  0.133  0.001
## +60             |  0.702  0.072  4.624 | -0.316  0.014 -2.078 | -0.033  0.000
##                 v.test  
## F               -1.930 |
## M                1.930 |
## employee         0.229 |
## middle           0.225 |
## non-worker      -0.091 |
## other worker     0.498 |
## senior           0.056 |
## student         -0.861 |
## workman          0.470 |
## +60             -0.218 |
## 
## Supplementary categorical variables (eta2)
##                   Dim.1 Dim.2 Dim.3  
## sex             | 0.014 0.000 0.012 |
## SPC             | 0.068 0.020 0.004 |
## age_Q           | 0.130 0.022 0.014 |
## frequency       | 0.011 0.007 0.026 |

5.2 Plotting and interpretation

Let’s start with a plot of the individuals (ind):

par(mfrow=c(2,2))
plot(mca, graph.type="ggplot", invisible=c("var","quali.sup","quanti.sup"),cex=0.8)

I see no real clusters, but there seem to be some duplicates. Dimension 1 covers 12.6% of the variance, and dimension 2 covers 12% of the variance, roughly the same. The scatter plot forms a triangle in this space, with a wide base forming around the Dim 1 = 0 and extending to the right.

Now let’s plot the var variables — Tea, price, and How — which the plot function for FactoMineR’s MCA colors black, red, and green, respectively.

plot(mca, graph.type="ggplot", invisible=c("ind","quali.sup","quanti.sup"), cex=0.8, habillage="quali")

Dimension 1 seems to cover: - most strongly (Tea: 0.61), whether the respondent prefers black tea or not (black +1.36, green -0.49, Earl Grey -0.44) - then (How: 0.334), whether the respondenr likes to add things to their tea (alone -0.28, lemon, milk, other all positive +) - and finally, no clear correlation with price, since both upscale and cheap are positive in dimension 1

Dimension 2 seems to cover: - whether the respondent is unlikely to prefer green tea (-1.636) - whether the respondent likes cheap tea (cheap +1.313, upscale -0.874)

The p_unknown category in the price has the highest coefficient here, which makes this one tough to interpret.

Finally, let’s look at the demographics of the individuals and see how they distribute over this space:

plot(mca, graph.type="ggplot", invisible=c("ind","var","quanti.sup"), cex=0.8, habillage="quali")

It seems that age is the best explaining variable for the differences in respondents tea drinking habits (dim 1: 0.130, dim 2: 0.022), followed by SPC, which does encode some of the same information as age, so perhaps it is redundant. If we interpret the dimensions from before we could roughly say that: - older people are more likely to prefer generic black tea to green tea or Earl Grey - younger people are more likely to prefer Earl Grey - green tea is most likely to be preferred by people in the 25-34 age range who are “employees” (office workers?) - older people prefer to add milk, lemon, or something else to their tea

There is also a slight difference in the distribution of answers between sexes: - men more likely to prefer black tea - women more likely to prefer Earl Grey or green tea

The frequency of tea drinking is all over the place and forms no clear trendline in the biplot.

1.3 Group mean and std. error plots

Now let’s do that same plot but instead of individuals, let’s plot the mean and standard error per group

std_rat_group_trends <-  rats_std %>% group_by(Time, Group) %>%
  summarise(mean = mean(std_weight), mean_gain = mean(std_weight_gain), se = sd(std_weight)/sqrt(n()),se_gain = sd(std_weight_gain)/sqrt(n())) %>%
  ungroup()

## `summarise()` has grouped output by 'Time'. You can override using the
## `.groups` argument.

sum_weight_plot <- std_rat_group_trends %>% ggplot(aes(x = Time, y = mean,ymin=mean-se, ymax=mean+se, color = Group)) +
  ggtitle("Bodyweight per group") +
  geom_line() +
  geom_point(size=3) +
  geom_linerange() +
  geom_ribbon(aes(fill=Group,alpha = 0.1)) +
  scale_y_continuous(name = "mean(std. weight) +/- se(std. weight)") +
  theme(legend.position="none")

sum_gain_plot <- std_rat_group_trends %>% ggplot(aes(x = Time,ymin=mean_gain-se_gain, ymax=mean_gain+se_gain, y = mean_gain, color = Group)) +
  ggtitle("Weight gain per group") +
  geom_line() +
  geom_point(size=3) +
  geom_linerange() +
  geom_ribbon(aes(fill=Group,alpha = 0.1)) +
  scale_y_continuous(name = "mean(std. weight gain) +/- se(std. weight gain)") +
  scale_alpha(guide="none") +
  theme(legend.position = c(0.95,0.08))

gridExtra::grid.arrange(sum_weight_plot,sum_gain_plot, ncol=2)

fig. 3.2, Bodyweight and weight gain over time per group

Looking at the weight gain plot specifically, we do see the three groups clearly separating, but we still don’t know how much of this effect would be because of the baseline weight of the rats or the treatment, since both of these variables differed between the groups.

1.4 Outliers

weight_boxplot <- final_rats %>% ggplot(aes(y=Bodyweight)) +
              geom_boxplot() +
              facet_wrap(~Group)

gain_boxplot <- final_rats %>% ggplot(aes(y=Weight_gain)) +
              geom_boxplot() +
              facet_wrap(~Group)

gridExtra::grid.arrange(weight_boxplot, gain_boxplot, ncol=2)

fig. 4.1, boxplots of the three groups

There is a rat in group 2 that seems to be increasing the variance in the group, which would make t-tests more difficult to run. Let’s remove it

outlier_ids <- (final_rats %>% filter(Bodyweight > 600))$ID
rats <- rats %>% filter(!(ID %in% outlier_ids))
initial_rats <- rats %>% filter(Time == initial_time)
final_rats <- rats %>% filter(Time == final_time)

weight_boxplot2 <- final_rats %>% ggplot(aes(y=Bodyweight)) +
              geom_boxplot() +
              facet_wrap(~Group)

gain_boxplot2 <- final_rats %>% ggplot(aes(y=Weight_gain)) +
              geom_boxplot() +
              facet_wrap(~Group)

gridExtra::grid.arrange(weight_boxplot2, gain_boxplot2, ncol=2)

fig. 4.2, boxplots of the three groups, outliers removed

There is still seemingly one outlier, but the difference in variances between the groups now seems better.

One issue is that the variance of group 1 is much smaller than groups 2 and 3 — likely because there are twice as many rats in that group.

1.5 Summary measure: final weight gain

Visually we have identified a between-groups effect, but there are two possible causes of the effect: the baseline weight and the treatment (diet).

Based on table 8.2 from MABS4IODS, I will use as summary measure the final value of the weight gain, because there is a growth trend in the data. Another option for summary measures would have been the regression coefficients of the bodyweight of each individual rat.

final_rats_12 <- final_rats %>% filter(Group != 3)
final_rats_13 <- final_rats %>% filter(Group != 2)
final_rats_23 <- final_rats %>% filter(Group != 1)

final_rats %>% group_by(Group) %>%
  summarise(var = var(Weight_gain))

## # A tibble: 3 × 2
##   Group    var
##   <fct>  <dbl>
## 1 1       86.1
## 2 2     1051  
## 3 3      326.

t.test(Weight_gain ~ Group, data = final_rats_12, var.equal = F)

## 
##  Welch Two Sample t-test
## 
## data:  Weight_gain by Group
## t = -2.0458, df = 2.1242, p-value = 0.1699
## alternative hypothesis: true difference in means between group 1 and group 2 is not equal to 0
## 95 percent confidence interval:
##  -116.22098   38.47098
## sample estimates:
## mean in group 1 mean in group 2 
##          23.125          62.000

t.test(Weight_gain ~ Group, data = final_rats_13, var.equal = F)

## 
##  Welch Two Sample t-test
## 
## data:  Weight_gain by Group
## t = -1.9138, df = 3.8172, p-value = 0.1316
## alternative hypothesis: true difference in means between group 1 and group 3 is not equal to 0
## 95 percent confidence interval:
##  -45.543111   8.793111
## sample estimates:
## mean in group 1 mean in group 3 
##          23.125          41.500

t.test(Weight_gain ~ Group, data = final_rats_23, var.equal = F)

## 
##  Welch Two Sample t-test
## 
## data:  Weight_gain by Group
## t = 0.98659, df = 2.932, p-value = 0.3981
## alternative hypothesis: true difference in means between group 2 and group 3 is not equal to 0
## 95 percent confidence interval:
##  -46.50239  87.50239
## sample estimates:
## mean in group 2 mean in group 3 
##            62.0            41.5

final_fit_gain <- lm(Weight_gain ~ Baseline_weight + Group, data = final_rats)
anova(final_fit_gain)

## Analysis of Variance Table
## 
## Response: Weight_gain
##                 Df Sum Sq Mean Sq F value   Pr(>F)   
## Baseline_weight  1 1308.3 1308.32  8.1039 0.015889 * 
## Group            2 4072.2 2036.11 12.6120 0.001423 **
## Residuals       11 1775.9  161.44                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

final_fit_gain_12 <- lm(Weight_gain ~ Baseline_weight + Group, data = final_rats_12)
anova(final_fit_gain_12)

## Analysis of Variance Table
## 
## Response: Weight_gain
##                 Df Sum Sq Mean Sq F value   Pr(>F)   
## Baseline_weight  1 2369.3 2369.30  15.279 0.004489 **
## Group            1 2392.3 2392.32  15.427 0.004370 **
## Residuals        8 1240.6  155.07                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

final_fit_gain_13 <- lm(Weight_gain ~ Baseline_weight + Group, data = final_rats_13)
anova(final_fit_gain_13)

## Analysis of Variance Table
## 
## Response: Weight_gain
##                 Df Sum Sq Mean Sq F value  Pr(>F)  
## Baseline_weight  1 662.01  662.01  6.9201 0.02733 *
## Group            1 957.27  957.27 10.0066 0.01149 *
## Residuals        9 860.98   95.66                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

final_fit_gain_23 <- lm(Weight_gain ~ Baseline_weight + Group, data = final_rats_23)
anova(final_fit_gain_23)

## Analysis of Variance Table
## 
## Response: Weight_gain
##                 Df Sum Sq Mean Sq F value  Pr(>F)  
## Baseline_weight  1 1870.8 1870.80  6.2693 0.06649 .
## Group            1  735.0  735.00  2.4631 0.19163  
## Residuals        4 1193.6  298.41                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

2.2 Independent linear model

Let’s first fit a model that ignores the correlation of repeated measures taken on the same subject.

bprs.lm <- lm(bprs ~ week + treatment, bprs)
summary(bprs.lm)

## 
## Call:
## lm(formula = bprs ~ week + treatment, data = bprs)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -22.454  -8.965  -3.196   7.002  50.244 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  46.4539     1.3670  33.982   <2e-16 ***
## week         -2.2704     0.2524  -8.995   <2e-16 ***
## treatment2    0.5722     1.3034   0.439    0.661    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 12.37 on 357 degrees of freedom
## Multiple R-squared:  0.1851, Adjusted R-squared:  0.1806 
## F-statistic: 40.55 on 2 and 357 DF,  p-value: < 2.2e-16

Based on the independent model, we find no statistically significant effect due to the treatment. This doesn’t necessarily mean anything, since the difference in intercepts and slopes between subjects is very high (as seen in figure 5.1) when compared to the difference in intercepts and slopes between groups.

The linear model does, however also identify the overall slope over time, which is -2.27.

2.3 Random intercept model

bprs.lmri <- lmer(bprs ~ week + treatment + (1|subject), data = bprs, REML = FALSE)
summary(bprs.lmri)

## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: bprs ~ week + treatment + (1 | subject)
##    Data: bprs
## 
##      AIC      BIC   logLik deviance df.resid 
##   2582.9   2602.3  -1286.5   2572.9      355 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.27506 -0.59909 -0.06104  0.44226  3.15835 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subject  (Intercept) 97.39    9.869   
##  Residual             54.23    7.364   
## Number of obs: 360, groups:  subject, 40
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  46.4539     2.3521  19.750
## week         -2.2704     0.1503 -15.104
## treatment2    0.5722     3.2159   0.178
## 
## Correlation of Fixed Effects:
##            (Intr) week  
## week       -0.256       
## treatment2 -0.684  0.000

The variance due to the random intercept is:

\(\hat{\sigma}^2_u = 97.38\)

And the residual variance is:

\(\hat{\sigma}^2 = 54.23\)

Since \(\hat{\sigma}^2_u \gt \hat{\sigma}^2\), I take this to mean that the random intercept was able to reduce the variance.

bprs$fitted_ri <- fitted(bprs.lmri)

ri_observed_plot <- ggplot(data=bprs, aes(x=week, color=treatment, group=subject )) +
    ggtitle("Observed values") +
    geom_line(aes(y=bprs)) +
    labs(color= "Treatment") +
    ylab("BPRS") +
    scale_fill_manual("", values =(c("se"="#cccccc"))) +
    scale_linetype_manual(guide="none") +
    scale_color_manual(values =(c("1"="#ff9999","2"="#9999ff"))) +
    scale_alpha(guide="none")

ri_fit_plot <- ggplot(data=bprs, aes(x=week, color=treatment, group=subject )) +
    ggtitle("Fitted values") +
    geom_line(aes(y=fitted_ri)) +
    labs(color= "Treatment") +
    ylab("BPRS") +
    scale_fill_manual("", values =(c("se"="#cccccc"))) +
    scale_linetype_manual(guide="none") +
    scale_color_manual(values =(c("1"="#ff9999","2"="#9999ff"))) +
    scale_alpha(guide="none")

gridExtra::grid.arrange(ri_observed_plot, ri_fit_plot, ncol=2)

fig. 6.1, Model with random intercept

confint(bprs.lmri)

## Computing profile confidence intervals ...

##                 2.5 %    97.5 %
## .sig01       7.918218 12.645375
## .sigma       6.828332  7.973573
## (Intercept) 41.744598 51.163179
## week        -2.565919 -1.974914
## treatment2  -5.885196  7.029641

2.3 Random intercept + slope model

bprs.lmris <- lmer(bprs ~ week + treatment + (week|subject), data = bprs, REML = FALSE)
summary(bprs.lmris)

## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: bprs ~ week + treatment + (week | subject)
##    Data: bprs
## 
##      AIC      BIC   logLik deviance df.resid 
##   2523.2   2550.4  -1254.6   2509.2      353 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4655 -0.5150 -0.0920  0.4347  3.7353 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  subject  (Intercept) 167.827  12.955        
##           week          2.331   1.527   -0.67
##  Residual              36.747   6.062        
## Number of obs: 360, groups:  subject, 40
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  45.9830     2.6470  17.372
## week         -2.2704     0.2713  -8.370
## treatment2    1.5139     3.1392   0.482
## 
## Correlation of Fixed Effects:
##            (Intr) week  
## week       -0.545       
## treatment2 -0.593  0.000

So now we have the following random effects:

\(\hat{\sigma}^2_u = 167.827\)

\(\hat{\sigma}^2_v = 2.331\)

And the residual variance is:

\(\hat{\sigma}^2 = 36.747\)

As the residual variance is less now, we already have a hint that the model is a better fit than the previous one.

Interestingly, the total variance is now much higher. Having discussed this with the course teachers, the intuitive reason for this is that the slopes being adjusted per-subject will tend to spread the intecepts more, as the lines can freely “rotate” into position now.

bprs$fitted_ris <- fitted(bprs.lmris)

ris_observed_plot <- ggplot(data=bprs, aes(x=week, color=treatment, group=subject )) +
    ggtitle("Observed values") +
    geom_line(aes(y=bprs)) +
    labs(color= "Treatment") +
    ylab("BPRS") +
    scale_fill_manual("", values =(c("se"="#cccccc"))) +
    scale_linetype_manual(guide="none") +
    scale_color_manual(values =(c("1"="#ff9999","2"="#9999ff"))) +
    scale_alpha(guide="none")

ris_fit_plot <- ggplot(data=bprs, aes(x=week, color=treatment, group=subject )) +
    ggtitle("Fitted values") +
    geom_line(aes(y=fitted_ris)) +
    labs(color= "Treatment") +
    ylab("BPRS") +
    scale_fill_manual("", values =(c("se"="#cccccc"))) +
    scale_linetype_manual(guide="none") +
    scale_color_manual(values =(c("1"="#ff9999","2"="#9999ff"))) +
    scale_alpha(guide="none")

gridExtra::grid.arrange(ris_observed_plot, ris_fit_plot, ncol=2)

fig. 6.2, Model with random intercept and slope

2.4 ANOVA between random intercept and random intercept + slope

anova(bprs.lmris,bprs.lmri)

## Data: bprs
## Models:
## bprs.lmri: bprs ~ week + treatment + (1 | subject)
## bprs.lmris: bprs ~ week + treatment + (week | subject)
##            npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## bprs.lmri     5 2582.9 2602.3 -1286.5   2572.9                         
## bprs.lmris    7 2523.2 2550.4 -1254.6   2509.2 63.663  2  1.499e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The deviance and AIC have both dropped, which tells us that adding the slope gives us a much better fit — AIC, the Akaike information criterion, is a measure of how much information is lost in the model, so less is better. At p=1.5e-14, we can confidently say that each subject definitely has a characteristic slope.

confint(bprs.lmris)

## Computing profile confidence intervals ...

##                  2.5 %     97.5 %
## .sig01      10.3102840 16.7099152
## .sig02      -0.8241817 -0.4128655
## .sig03       1.1562015  2.0277746
## .sigma       5.5925524  6.6009205
## (Intercept) 40.5849755 51.2468100
## week        -2.8150999 -1.7257334
## treatment2  -4.9313699  7.9592153

2.5 Random intercept + slope with interaction

bprs.lmris_int <- lmer(bprs ~ week + treatment + week * treatment + (week|subject) , data = bprs, REML = FALSE)
summary(bprs.lmris_int)

## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: bprs ~ week + treatment + week * treatment + (week | subject)
##    Data: bprs
## 
##      AIC      BIC   logLik deviance df.resid 
##   2523.5   2554.5  -1253.7   2507.5      352 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4747 -0.5256 -0.0866  0.4435  3.7884 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  subject  (Intercept) 164.204  12.814        
##           week          2.203   1.484   -0.66
##  Residual              36.748   6.062        
## Number of obs: 360, groups:  subject, 40
## 
## Fixed effects:
##                 Estimate Std. Error t value
## (Intercept)      47.8856     2.9840  16.047
## week             -2.6283     0.3752  -7.006
## treatment2       -2.2911     4.2200  -0.543
## week:treatment2   0.7158     0.5306   1.349
## 
## Correlation of Fixed Effects:
##             (Intr) week   trtmn2
## week        -0.668              
## treatment2  -0.707  0.473       
## wek:trtmnt2  0.473 -0.707 -0.668

The variances have barely changed, and the plot below also looks almost exactly the same as the previous model.

bprs$fitted_ris_int <- fitted(bprs.lmris_int)

int_observed_plot <- ggplot(data=bprs, aes(x=week, color=treatment, group=subject )) +
    ggtitle("Observed values") +
    geom_line(aes(y=bprs)) +
    labs(color= "Treatment") +
    ylab("BPRS") +
    scale_fill_manual("", values =(c("se"="#cccccc"))) +
    scale_linetype_manual(guide="none") +
    scale_color_manual(values =(c("1"="#ff9999","2"="#9999ff"))) +
    scale_alpha(guide="none")

int_fit_plot <- ggplot(data=bprs, aes(x=week, color=treatment, group=subject )) +
    ggtitle("Fitted values") +
    geom_line(aes(y=fitted_ris_int)) +
    labs(color= "Treatment") +
    ylab("BPRS") +
    scale_fill_manual("", values =(c("se"="#cccccc"))) +
    scale_linetype_manual(guide="none") +
    scale_color_manual(values =(c("1"="#ff9999","2"="#9999ff"))) +
    scale_alpha(guide="none")

gridExtra::grid.arrange(int_observed_plot, int_fit_plot, ncol=2)

fig. 6.3, Model with random intercept and slope, including interaction between week and treatment

2.6 ANOVA between rand. intercept + slope model with interaction and without

anova(bprs.lmris_int,bprs.lmris)

## Data: bprs
## Models:
## bprs.lmris: bprs ~ week + treatment + (week | subject)
## bprs.lmris_int: bprs ~ week + treatment + week * treatment + (week | subject)
##                npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)
## bprs.lmris        7 2523.2 2550.4 -1254.6   2509.2                    
## bprs.lmris_int    8 2523.5 2554.6 -1253.7   2507.5  1.78  1     0.1821

We see that there is a negligible drop in deviance and even an increase in AIC, and p=0.18 — there is not a significant effect due to a week-treatment interaction.

confint(bprs.lmris_int)

## Computing profile confidence intervals ...

##                       2.5 %     97.5 %
## .sig01           10.2191215 16.4882552
## .sig02           -0.8181446 -0.4021099
## .sig03            1.1186282  1.9764648
## .sigma            5.5925518  6.6009202
## (Intercept)      41.8936922 53.8774190
## week             -3.3816817 -1.8749850
## treatment2      -10.7648856  6.1826634
## week:treatment2  -0.3495621  1.7812288

2.7 Interpretation

We have to conclude that, using a linear mixed effects model, we cannot find a significant effect due to the treatment.